#include "ParticleSystem.h"
#include "ParticleSystem.cuh"

#include <thrust/scan.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/device_ptr.h>
#include <thrust/sort.h>

#include <stdlib.h>
#include <stdio.h>

#define NTHREADS 512

texture<float> texPos;
texture<float> texVel;
texture<float> texForce;

__device__ __constant__ float devRadius[NUM_PARTICLES];

/******************************************************************************
    AUXILIARY FUNCTIONS FORWARD DECLARATIONS
******************************************************************************/
void buildGrid( float3* devPosition, int numParticles,
                thrust::device_vector<int>& pairPart,
                thrust::device_vector<int>& cellLimitBegin,
                thrust::device_vector<int>& cellLimitEnd   );

__host__ void handleError( const char* filename, const int lineNumber, cudaError_t error )
{
    if (error != cudaSuccess)
    {
        printf("%s:%i: ", filename, lineNumber);
    }
        
    switch (error)
    {
    case cudaSuccess: // do nothing
        break;
   
    case cudaErrorMissingConfiguration:
        printf("cudaErrorMissingConfiguration: ");
        break;    
        
    case cudaErrorMemoryAllocation:
        printf("cudaErrorMemoryAllocation: ");
        break;
            
    case cudaErrorInitializationError:
        printf("cudaErrorInitializationError: ");
        break;
            
    case cudaErrorLaunchFailure:
        printf("cudaErrorLaunchFailure: ");
        break;
    
    case cudaErrorPriorLaunchFailure:
        printf("cudaErrorPriorLaunchFailure: ");
        break;
        
    case cudaErrorLaunchTimeout:
        printf("cudaErrorLaunchTimeout: ");
        break;
        
    case cudaErrorLaunchOutOfResources:
        printf("cudaErrorLaunchOutOfResources: ");
        break;
        
    case cudaErrorInvalidDeviceFunction:
        printf("cudaErrorInvalidDeviceFunction: ");
        break;
        
    case cudaErrorInvalidConfiguration:
        printf("cudaErrorInvalidConfiguration: ");
        break;
        
    case cudaErrorInvalidDevice:
        printf("cudaErrorInvalidDevice: ");
        break;
        
    case cudaErrorInvalidValue:
        printf("cudaErrorInvalidValue: ");
        break;
        
    case cudaErrorInvalidPitchValue:
        printf("cudaErrorInvalidPitchValue: ");
        break;
        
    case cudaErrorInvalidSymbol:
        printf("cudaErrorInvalidSymbol: ");
        break;
    
    case cudaErrorMapBufferObjectFailed:
        printf("cudaErrorMapBufferObjectFailed: ");
        break;
        
    case cudaErrorUnmapBufferObjectFailed:
        printf("cudaErrorUnmapBufferObjectFailed: ");
        break;
        
    case cudaErrorInvalidHostPointer:
        printf("cudaErrorInvalidHostPointer: ");
        break;
        
    case cudaErrorInvalidDevicePointer:
        printf("cudaErrorInvalidDevicePointer: ");
        break;
        
    case cudaErrorInvalidTexture:
        printf("cudaErrorInvalidTexture: ");
        break;
        
    case cudaErrorInvalidTextureBinding:
        printf("cudaErrorInvalidTextureBinding: ");
        break;
        
    case cudaErrorInvalidChannelDescriptor:
        printf("cudaErrorInvalidChannelDescriptor: ");
        break;
        
    case cudaErrorInvalidMemcpyDirection:
        printf("cudaErrorInvalidMemcpyDirection: ");
        break;
        
    case cudaErrorAddressOfConstant:
        printf("cudaErrorAddressOfConstant: ");
        break;
        
    case cudaErrorTextureFetchFailed:
        printf("cudaErrorTextureFetchFailed: ");
        break;
        
    case cudaErrorTextureNotBound:
        printf("cudaErrorTextureNotBound: ");
        break;
        
    case cudaErrorSynchronizationError:
        printf("cudaErrorSynchronizationError: ");
        break;
        
    case cudaErrorInvalidFilterSetting:
        printf("cudaErrorInvalidFilterSetting: ");
        break;
        
    case cudaErrorInvalidNormSetting:
        printf("cudaErrorInvalidNormSetting: ");
        break;
        
    case cudaErrorMixedDeviceExecution:
        printf("cudaErrorMixedDeviceExecution: ");
        break;
        
    case cudaErrorCudartUnloading:
        printf("cudaErrorCudartUnloading: ");
        break;
        
    case cudaErrorUnknown:
        printf("cudaErrorUnknown: ");
        break;
        
    case cudaErrorNotYetImplemented:
        printf("cudaErrorNotYetImplemented: ");
        break;
        
    case cudaErrorMemoryValueTooLarge:
        printf("cudaErrorMemoryValueTooLarge: ");
        break;
        
    case cudaErrorInvalidResourceHandle:
        printf("cudaErrorInvalidResourceHandle: ");
        break;
        
    case cudaErrorNotReady:
        printf("cudaErrorNotReady: ");
        break;
        
    case cudaErrorInsufficientDriver:
        printf("cudaErrorInsufficientDriver: ");
        break;
        
    case cudaErrorSetOnActiveProcess:
        printf("cudaErrorSetOnActiveProcess: ");
        break;
        
    case cudaErrorNoDevice:
        printf("cudaErrorNoDevice: ");
        break;
        
    case cudaErrorStartupFailure:
        printf("cudaErrorStartupFailure: ");
        break;
        
    case cudaErrorApiFailureBase:
        printf("cudaErrorApiFailureBase: ");
        break;
    }
    
  	if(error != cudaSuccess)
  	{
   	 // print the CUDA error message and exit
  	  printf("%s\n", cudaGetErrorString(error));
  	  exit(-1);
  	}
}

/******************************************************************************
    DEVICE KERNELS
******************************************************************************/
__device__ int makeCellId( int xcell, int ycell, int zcell )
{
    return ycell*CELL_NUM*CELL_NUM + zcell*CELL_NUM + xcell;
}
__device__ int makeCellId( float3 pos )
{
    // float position in the grid.
    pos = ( pos - BOX_ORIGIN_VEC ) / CELL_DIM;

    // casting to int will round down.
    // find the cell id of the rounded down grid coords
    return makeCellId( (int)pos.x, (int)pos.y, (int)pos.z );
}
__device__ float3 getCollisionVector( float3 posA, float radiusA, float3 posB, float radiusB )
{
    //calculates vector from A to B
    float3 collisionVec = posB - posA;

    float l = length(collisionVec);
    //if they are interpenetrating
    if( l < (radiusA + radiusB) )
    {
        return ((radiusA+radiusB)*normalize(collisionVec)) - collisionVec;
    }
    
    return make_float3( 0.0f, 0.0f, 0.0f );
}

__device__ float getRelativeVelocity( float3 velA, float3 velB, float3 dir_AB )
{
    return dot( (velA-velB), dir_AB );
}


__device__ void getOverlappingCells( float3* devPos, int pid, int* cellIds )
{
    // get position of the current particle
    float3 pos = make_float3(tex1Dfetch( texPos  , 3*pid   ),
                             tex1Dfetch( texPos  , 3*pid+1 ),
                             tex1Dfetch( texPos  , 3*pid+2 ));
    
    // get maximum and minimum x,y,z possible values in the grid
    float3 posMax = (pos + devRadius[pid] - BOX_ORIGIN_VEC)/CELL_DIM;
    float3 posMin = (pos - devRadius[pid] - BOX_ORIGIN_VEC)/CELL_DIM;

    // get maximum and minimum xyz of the cells
    int3 cellMax = make_int3( posMax.x, posMax.y, posMax.z );
    int3 cellMin = make_int3( posMin.x, posMin.y, posMin.z );

    // run through all the possible xyz combinations, and add the cell ids to the array
    int i = 0;
    for( int xcell = cellMin.x; xcell <= cellMax.x; xcell++ )
    {
        for( int ycell = cellMin.y; ycell <= cellMax.y; ycell++ )
        {
            for( int zcell = cellMin.z; zcell <= cellMax.z; zcell++ )
            {
                cellIds[i++] = makeCellId( xcell, ycell, zcell );
            }
        }
    }
}

__device__ float3 getCollisionPoint( float3 posA, float radiusA, float3 posB, float radiusB )
{
    // P1 + |d|*(r1-0.5*(r1+r2-d))
    return posA + normalize(posB-posA)*( radiusA - 0.5f * (radiusA+radiusB-length(posB-posA)) );
}
/******************************************************************************
    GLOBAL KERNELS
******************************************************************************/
__global__ void kernelCountPairs( float3* devPos, int numParticles, int* devPairCountArray )
{
    // get particle id from thread info
    int pid = threadIdx.x + blockIdx.x * blockDim.x;
    if( pid >= numParticles ) return;

    // fetch position from texture
    float3 position = make_float3(tex1Dfetch( texPos  , 3*pid   ),
                                  tex1Dfetch( texPos  , 3*pid+1 ),
                                  tex1Dfetch( texPos  , 3*pid+2 ));
    
    // calculate max and minimum grid indexes
    int3 maxCell = make_int3((position.x + devRadius[pid] - BOX_ORIGIN_X)/CELL_DIM,
                             (position.y + devRadius[pid] - BOX_ORIGIN_Y)/CELL_DIM,
                             (position.z + devRadius[pid] - BOX_ORIGIN_Z)/CELL_DIM);

    int3 minCell = make_int3((position.x - devRadius[pid] - BOX_ORIGIN_X)/CELL_DIM,
                             (position.y - devRadius[pid] - BOX_ORIGIN_Y)/CELL_DIM,
                             (position.z - devRadius[pid] - BOX_ORIGIN_Z)/CELL_DIM);

    // number of cells overlapped in 1 dimension -> max-min+1
    // total number of pair equal -> n_overlaps.x * n_overlaps.y * n_overlaps.z
    devPairCountArray[pid] = (maxCell.x - minCell.x + 1) * (maxCell.y - minCell.y + 1) * (maxCell.z - minCell.z + 1);
}

__global__ void kernelFillPairArrays( float3* devPos, int numParticles, int* devCountScan, int* pairPart, int* pairCell )
{
    // get the particle id from thread info
    int pid = threadIdx.x + blockIdx.x * blockDim.x;
    if( pid >= numParticles ) return;

    // create a vector to hold the ids of the cells this particle overlaps
    int overlapCells[MAX_OVERLAPS];
    getOverlappingCells( devPos, pid, overlapCells );

    // find the indexes of the pair vectors in wich I need to write on
    int begin = devCountScan[pid];
    int end   = devCountScan[pid+1];

    // iterate filling the pair vectors
    for( int i=begin; i<end; ++i )
    {
        pairPart[i] = pid;
        pairCell[i] = overlapCells[i-begin];
    }
}

__global__ void kernelFillCellLimitArrays( int* pairCell ,int nPairs, int* cellLimitBegin , int* cellLimitEnd )
{
    int pairId = threadIdx.x + blockIdx.x * blockDim.x;

    if( pairId >= nPairs ) return;

    int cellId = pairCell[pairId];

    //if( pairId >= nPairs ) return;

    if( pairId == 0 ) //Treat first case
    {
        for( int i = 0; i < cellId; ++i )
        {
            cellLimitBegin[i] = 0;
            cellLimitEnd[i] = 0;
        }

        cellLimitBegin[cellId] = 0;

        return;
    }

    if( pairCell[pairId] > pairCell[pairId-1] )
    {
        
        for( int i = pairCell[pairId-1]; i<cellId; ++i )
        {
            cellLimitEnd[i] = pairId;
            cellLimitBegin[i+1] = pairId;
        }
        //cellLimitEnd[cellId-1] = pairId;
        ///cellLimitBegin[cellId] = pairId;
    }

    if( (pairId+1) == nPairs ) //Treat last case
    {
        cellLimitEnd[cellId] = nPairs;

        for( int i=cellId+1; i<CELL_TOTAL; ++i )
        {
            cellLimitBegin[i] = nPairs;
            cellLimitEnd[i] = nPairs;
        }
    }
}

__global__ void kernelCollideParticles( float dt, int numParticles,
                                        float3* devPos, float3* devVel, float3* devForce,
                                        float3* outPos, float3* outVel, float3* outForce,
                                        int* particlesOrdered, int* cellLimitBegin, int* cellLimitEnd )
{
    int pid = threadIdx.x + blockIdx.x * blockDim.x;
    if( pid >= numParticles ) return;

    float3 devPosPid   = make_float3(tex1Dfetch( texPos  , 3*pid   ),
                                     tex1Dfetch( texPos  , 3*pid+1 ),
                                     tex1Dfetch( texPos  , 3*pid+2 ));
    float3 devVelPid   = make_float3(tex1Dfetch( texVel  , 3*pid   ),
                                     tex1Dfetch( texVel  , 3*pid+1 ),
                                     tex1Dfetch( texVel  , 3*pid+2 ));
    float3 devForcePid = make_float3(tex1Dfetch( texForce, 3*pid   ),
                                     tex1Dfetch( texForce, 3*pid+1 ),
                                     tex1Dfetch( texForce, 3*pid+2 ));
    float radiusPid    = devRadius[pid];
    float3 inPos = devPosPid;
    float3 inVel = devVelPid;


    //add force of gravity
    devForcePid.y -= GRAVITY * 1;

    //Start collision Treatment
    int vecCellIDs[MAX_OVERLAPS];

    for( int i=0; i<MAX_OVERLAPS; ++i )
    {
        vecCellIDs[i] = -1;
    }

    getOverlappingCells( devPos, pid, vecCellIDs );

    // go into each overlapping cell to check for collisions
    for( int c = 0; c < MAX_OVERLAPS; ++c )
    {
        if( vecCellIDs[c] < 0 ) break;

        int cellBegin = cellLimitBegin[vecCellIDs[c]];
        int cellEnd   = cellLimitEnd[vecCellIDs[c]];

        // check collision with each particle of the cell
        for( int i = cellBegin; i < cellEnd; ++i )
        {
            int pidOther = particlesOrdered[i];

            //do not collide it with itself!
            if( pidOther == pid ) continue;

            float3 devPosPidOther = make_float3(tex1Dfetch( texPos, 3*pidOther   ),
                                                tex1Dfetch( texPos, 3*pidOther+1 ),
                                                tex1Dfetch( texPos, 3*pidOther+2 ));
            float3 devVelPidOther = make_float3(tex1Dfetch( texVel, 3*pidOther   ),
                                                tex1Dfetch( texVel, 3*pidOther+1 ),
                                                tex1Dfetch( texVel, 3*pidOther+2 ));

            //get the collision vector between particles.
            float3 collisionVec = getCollisionVector(inPos, radiusPid, devPosPidOther, devRadius[pidOther]);

            //if it is a null vector, fetch another pair of particles.
            if( length(collisionVec) < 0.001f ) continue;

            //if the collision point did not happen inside this cell, fetch another pair of particles.
            float3 collisionPoint = getCollisionPoint(inPos, radiusPid, devPosPidOther, devRadius[pidOther]);
            if( makeCellId(collisionPoint) != vecCellIDs[c] ) continue;

            //change the positions to a valid state
            devPosPid -= collisionVec/2;

            //get the direction of the collision
            float3 collisionDir = normalize(collisionVec);

            //find relative velocity between particles.
            float velRel = getRelativeVelocity(inVel, devVelPidOther, collisionDir);

            //calculate the impulse factor ( j )
            //j = (-(1+EPSILON)*vrel) / ((1/ma)+(1/mb))
            //rearranging so we have less divisions
            //j = 0.5 * (-1-EPSILON) * vrel * ma * mb
            //calculate the impulse vector ( J )
            //J = j * direction of the collision
            float3 J = collisionDir * (0.5f * (-1-RESTITUTION) * velRel * 1 * 1);

            //add the impulse to the particles
            //va2 - va1 = J/ma   and    vb2 - vb1 = -J/mb
            //rearranging
            //va2 = va1 + J/ma   and    vb2 = vb1 - J/mb
            devVelPid += J/1;
        }
    }

    //Treat collision with the floor
    if( devPosPid.y < radiusPid )
    {
        devPosPid.y = radiusPid;
        devVelPid.y *= -RESTITUTION;
    }

    //Treat collision with the ceilling
    if( devPosPid.y > BOX_DIM-radiusPid )
    {
        devPosPid.y = BOX_DIM-radiusPid;
        devVelPid.y *= -RESTITUTION;
    }

    //Treat collision with the walls
    if( devPosPid.x > (BOX_DIM/2-radiusPid) )
    {
        devPosPid.x = BOX_DIM/2-radiusPid;
        devVelPid.x *= -RESTITUTION;
    }
    if( devPosPid.x < -(BOX_DIM/2-radiusPid) )
    {
        devPosPid.x = -(BOX_DIM/2-radiusPid);
        devVelPid.x *= -RESTITUTION;
    }
    if( devPosPid.z < -(BOX_DIM/2-radiusPid) )
    {
        devPosPid.z = -(BOX_DIM/2-radiusPid);
        devVelPid.z *= -RESTITUTION;
    }
    if( devPosPid.z > BOX_DIM/2-radiusPid )
    {
        devPosPid.z = BOX_DIM/2-radiusPid;
        devVelPid.z *= -RESTITUTION;
    }

    // write to Global Memory
    outPos[pid]   = devPosPid;
    outVel[pid]   = devVelPid;
    outForce[pid] = devForcePid;
}

__global__ void kernelUpdateEuler( float dt, int numParticles, float3* devPos, float3* devVel, float3* devForce )
{
    int pid = threadIdx.x + blockIdx.x * blockDim.x;
    if( pid >= numParticles ) return;

    float3 nowPos = devPos[pid];
    float3 nowVel = devVel[pid];

    //Update Euler
    devPos[pid] = nowPos + (nowVel * dt);
    devVel[pid] = nowVel + (devForce[pid] * dt);

    //clean forces
    devForce[pid] = make_float3(0.0f);
}

/******************************************************************************
    CLASS METHODS DEFINITION
******************************************************************************/
void ParticleSystem::updateSpacePartitionCuda( float dt )
{
    static thrust::device_vector<int> particlesOrdered( NUM_PARTICLES * MAX_OVERLAPS );
    static thrust::device_vector<int> cellLimitBegin(CELL_TOTAL);
    static thrust::device_vector<int> cellLimitEnd(CELL_TOTAL);

    float3 *inPos , *inVel , *inForce;
    float3 *outPos, *outVel, *outForce;

    if( buffOneInput )
    {
        inPos    = devPosition;
        inVel    = devVelocity;
        inForce  = devForce;

        outPos   = devPosition_2;
        outVel   = devVelocity_2;
        outForce = devForce_2;
    }
    else
    {
        inPos    = devPosition_2;
        inVel    = devVelocity_2;
        inForce  = devForce_2;

        outPos   = devPosition;
        outVel   = devVelocity;
        outForce = devForce;
    }

    cudaBindTexture( NULL, texPos  , inPos  , sizeof( float ) * NUM_PARTICLES * 3 );
    cudaBindTexture( NULL, texVel  , inVel  , sizeof( float ) * NUM_PARTICLES * 3 );
    cudaBindTexture( NULL, texForce, inForce, sizeof( float ) * NUM_PARTICLES * 3 );


    kernelUpdateEuler<<<(numParticles+NTHREADS-1)/NTHREADS, NTHREADS>>>( dt, numParticles, inPos, inVel, inForce );

    buildGrid( inPos, numParticles, particlesOrdered, cellLimitBegin, cellLimitEnd );

    kernelCollideParticles<<<(numParticles+NTHREADS-1)/NTHREADS, NTHREADS>>>(dt, numParticles, inPos, inVel, inForce,
                                                                             outPos, outVel, outForce,
                                                                             thrust::raw_pointer_cast(&particlesOrdered[0]),
                                                                             thrust::raw_pointer_cast(&cellLimitBegin[0]),
                                                                             thrust::raw_pointer_cast(&cellLimitEnd[0])      );
    buffOneInput = !(buffOneInput);
}

void ParticleSystem::initCudaArrays()
{
    // init position arrays
    cudaMalloc( (void**)&devPosition  , sizeof( float3 ) * NUM_PARTICLES );
    cudaMalloc( (void**)&devPosition_2, sizeof( float3 ) * NUM_PARTICLES );
    cudaMemcpy( devPosition  , pPosition, sizeof( float3 ) * NUM_PARTICLES, cudaMemcpyHostToDevice );
    cudaMemcpy( devPosition_2, pPosition, sizeof( float3 ) * NUM_PARTICLES, cudaMemcpyHostToDevice );

    // init velocity arrays
    cudaMalloc( (void**)&devVelocity  , sizeof( float3 ) * NUM_PARTICLES );
    cudaMalloc( (void**)&devVelocity_2, sizeof( float3 ) * NUM_PARTICLES );
    cudaMemcpy( devVelocity  , pVelocity, sizeof( float3 ) * NUM_PARTICLES, cudaMemcpyHostToDevice);
    cudaMemcpy( devVelocity_2, pVelocity, sizeof( float3 ) * NUM_PARTICLES, cudaMemcpyHostToDevice);

    // init force arrays
    cudaMalloc( (void**)&devForce  , sizeof( float3 ) * NUM_PARTICLES );
    cudaMalloc( (void**)&devForce_2, sizeof( float3 ) * NUM_PARTICLES );
    cudaMemcpy( devForce  , pForce, sizeof( float3 ) * NUM_PARTICLES, cudaMemcpyHostToDevice );
    cudaMemcpy( devForce_2, pForce, sizeof( float3 ) * NUM_PARTICLES, cudaMemcpyHostToDevice );

    // init radius array
    cudaMemcpyToSymbol( "devRadius", pRadius, sizeof(float) * NUM_PARTICLES );
}

void ParticleSystem::update(float dt)
{
    updateSpacePartitionCuda( dt );
#ifdef RENDER
    cudaMemcpy( pPosition,  devPosition, sizeof(float3)*numParticles, cudaMemcpyDeviceToHost );
#endif
}

/******************************************************************************
    AUXILIARY FUNCTIONS
******************************************************************************/
void buildGrid( float3* devPosition, int numParticles,
                thrust::device_vector<int>& pairPart,
                thrust::device_vector<int>& cellLimitBegin,
                thrust::device_vector<int>& cellLimitEnd   )
{
    //Calculate the number of pairs per particle (particle, cell)
    static thrust::device_vector<int> devCount(NUM_PARTICLES);
    kernelCountPairs<<<(numParticles+NTHREADS-1)/NTHREADS, NTHREADS>>>(devPosition, numParticles, thrust::raw_pointer_cast(&devCount[0]) );

    #ifdef PRINT_DEBUG
        thrust::host_vector<int> host_devCount = devCount;
        for( int i=0; i<NUM_PARTICLES; ++i )
        {
            printf( "particle %d overlaps %d cells\n", i, host_devCount[i] );
        }
    #endif
    
    //Exclusive scan (in place)
    static thrust::device_vector<int> devCountScan(NUM_PARTICLES+1);
    thrust::exclusive_scan(devCount.begin(), devCount.end(), devCountScan.begin());
    
    //Find total number of pairs
    int totalPairs = devCount[numParticles-1] + devCountScan[numParticles-1]; 
    devCountScan[numParticles] = totalPairs;

    #ifdef PRINT_DEBUG
        thrust::host_vector<int> host_devCountScan = devCountScan;
        for( int i=0; i<=NUM_PARTICLES; ++i )
        {
            printf( "devCountScan[%d] : %d\n", i, host_devCountScan[i] );
        }
    #endif

    //Create vector of pairs, which are in fact, two vectors.
    thrust::device_vector<int> pairCell( totalPairs );
    kernelFillPairArrays<<<(numParticles+NTHREADS-1)/NTHREADS, NTHREADS>>>( devPosition, numParticles,
                                                thrust::raw_pointer_cast(&devCountScan[0]),
                                                thrust::raw_pointer_cast(&pairPart[0]),
                                                thrust::raw_pointer_cast(&pairCell[0]) );
        
    #ifdef PRINT_DEBUG
        thrust::host_vector<int> host_pairPart = pairPart;
        thrust::host_vector<int> host_pairCell = pairCell;
        for( int i=0; i<totalPairs; ++i )
        {
            printf( "Particle %d is in cell %d\n", host_pairPart[i], host_pairCell[i] );
        }
    #endif

    //Now sort the vectors by cell index
    thrust::sort_by_key(pairCell.begin(), pairCell.end(), pairPart.begin());

    #ifdef PRINT_DEBUG
        thrust::host_vector<int> host_pairPart2 = pairPart;
        thrust::host_vector<int> host_pairCell2 = pairCell;
        for( int i=0; i<totalPairs; ++i )
        {
            printf( "%d. Cell %d contains particle %d\n", i, host_pairCell2[i], host_pairPart2[i] );
        }
    #endif

    //Generate the vector of cell limits. They indicate the index where the cell begins and ends in the pair vector
    kernelFillCellLimitArrays<<<(totalPairs+NTHREADS-1)/NTHREADS, NTHREADS>>>(thrust::raw_pointer_cast(&pairCell[0]),
                                                                              totalPairs,
                                                                              thrust::raw_pointer_cast(&cellLimitBegin[0]),
                                                                              thrust::raw_pointer_cast(&cellLimitEnd[0]) );

    #ifdef PRINT_DEBUG
        thrust::host_vector<int> auxVecBegin = cellLimitBegin;
        thrust::host_vector<int> auxVecEnd = cellLimitEnd;
        printf( "total pairs: %d\n", totalPairs );
        for( int i=0; i<CELL_TOTAL; ++i )
        {
            printf("cell %d: (%d,%d)\n", i, auxVecBegin[i], auxVecEnd[i]);
        }
    #endif
}